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We report extensive self-consistent calculations of jellium surface energies, by going beyond the local- 
density approximation. The density-response function of a bounded free-electron gas is evaluated 
within the random-phase approximation, with use of self-consistent electron density profiles. The 
exchange- correlation energy is then obtained from an exact adiabatic fluctuation-dissipation formula. 
We also investigate quantum-size effects and the extrapolation of finite-slab calculations to the 
infinite- width limit. 



I. INTRODUCTION 

Density-functional theory (DFT),§ in the local-density 
approximation (LDA) and its various gradient-corrected 
forms,a has been the most vifidely used |£aethod to inves- 
tigate the properties of solid surfaces.0Ll However, the 
rapidly varying electron density near the surface does 
not justify a priori the use of local or quasilocal approx- 
imations. In particular, these approximations are well 
known to yield an inaccurate description of the position- 
dependent exchange-correlation (xc) hole density, they 
lead to a surface barrier which does not present the cor- 
rect image-like asymptotic behavior, and there has been a 
long-standing controversy about the extent to wiijct the 
LDA accounts for the exact xc surface energy.l3iid Re- 
cently, the wave-function-based Fermi hypernetted-chain 
(FHNC)E3 and quantum Monte Carlo (QMC)^^ predic- 
tions have been found to disagree with modern density- 
functional calculations of jellium surface, energies that go 
beyond a local-density approximation.EjtS 

In this paper, we present extensive self-consistent cal- 
culations of jellium surface energies. The dynamical 
density-response function of a bounded free electron gas 
(FEG) is evaluated within the random-phase approxima- 
tion (RPA) Jla from the knowledge of the non- interacting 
density-response function. The xc energy of jellium slabs 
is then obtained via a coupling-constant integration over 
the interacting density-response function. For compari- 
son, we first coirsider a simple non-sclfconsistent micro- 
scopic model of the surface, the so-called infinite-barrier 
model (IBM), for which analytical insight is possible 
by virtue of the simple nature of the one-electron wave 
functions. Full self-consistent calculations are performed 
by employing either Kohn-Sham or Hartree orbitals to 
construct the non-interacting density-response function. 
LDA calculations of jellium surface energies are also car- 
ried out, within the same density-response framework, 
and a comparison between our local and non-local calcu- 



lations indicate that the error introduced by the LDA is 
small. r— in 

Quantum-size effects (QSE)ElEl on the surface energy 
of jellium slabs are also investigated, and the extrapola- 
tion of finite-slab calculations to the infinite-width limit 
is discussed. This issue is important, in view of the po- 
tential application of this procedure to the case of more 
complex slab calculations involving the use of the band 
structure and/or QMC simulations. 



II. THE JELLIUM SURFACE ENERGY 

We consider a jellium slab of thickness a normal to the 
z axis, consisting of a fixed uniform positive background 
of density 



n+{z) 



n, -a/2<z<a/2 
0, elsewhere, 



(1) 



plus a neutralizing cloud of interacting electrons. 

The surface energy is defined as the energy per unit 
area required to split the solid into two separate halves 
along a plane. For a solid that is translationally invariant 
in the plane of the surface, DFT shows thatQ 



a = CTs + (Jes + CTxc, 



(2) 



where as represents the kinetic surface energy of non- 
interacting Kohn-Sham electrons 
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dz n{z) {veff [n] (z) - v^ff [n] (0)} , (3) 



aes is the electrostatic surface energy of all positive and 
negative charge distributions 



1 



dz [n{z) — nj^{z)] ip{z), 



and a^c is the xc surface energy 

dzn{z) {£,,[n](z)-£™/(n)} 



Here, n{z) is the electron density 
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Vef f [n] (z) is the Kohn-Sham effective one-electron poten- 
tial 



VeffMiz) = ip{z) + Vxc[n]{z), 



(7) 



f{z) and Ua;c['T-](2) being the electrostatic and xc poten- 
tial, respectively, ea;cM(2:) represents an xc energy per 
particle at point z, and e^™^ {n) is the xc energy per par- 
ticle of a uniform electron gas of density n. Si and 4'i{z) 
are eigenvalues and eigenfunctions of the one-dimensional 
Kohn-Sham hamiltonian describing motion normal to the 
surface, Im denotes the highest occupied level, Ep rep- 
resents the Fermi level, and qp = (37r^n)^/'^ is the Fermi 
momentum. In the LDA, the xc surface energy is ob- 
tained by simply replacing ea;c[n](z) by the xc energy per 
particle of a uniform electron gas of density n(z), i.e., 
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dzn{z) {eT:f[n{z)]-elT^{n)]. (8) 



The exact xc energy per particle at point r can be 
obtained as the energy due to the interaction between an 
electron at r and its averaged xc hole: 



exc[n][v) 



dr 
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(9) 



where fixd'^TY') is defined by adiabatically switching on 
the electron-electron interaction via a coupling constant 
A, i.e., v^{v — r') = Ae^/|r — r'| and by adding, at the 
same time, an external potential so as to maintain the 
true (A — 1) ground-state density in the presencCj-pf the 
modified electron-electron interaction. One writesO 



n^c(r, r') = / dAn^c(r,r'), 



(10) 



where ri^ci''^, r') is the xc-hole density of a fictitious 
system at coupling jStrength A. Using the fluctuation- 
dissipation theorem,E3 one may write 



=(r,r') = 



1 n 

-n(r) (5(r - r')] 



dujx'{v,v']iuj) 
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where x'*'(i"; i"'; ^) represents the so-called density- 
response function of the electron system. 



Time-dependent DFT (TDDFT) shows that the exoft 
density-response function satisfies the integral equationE^ 



x\r,r';w) =X°(i-,i-';^)+ j dri J dra 

xx°(r, ri;w){«^(ri,r2)+/i',[n](ri,r2;^)} 
xx^(r2,r';w). 



where 
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v^^[n]{r,uA. being the exact frequency-dependent xc 
potential,Efl and x°(r, r'; ij^), the density-response func- 
tion of non- interacting Kohn-Sham electrons, i.e., in- 
dependent electrons moving in the effective potential 
Veff[n]{r) entering the Kohn-Sham equation of DFT. 

If the interacting density-response function x'^(r, r'; w) 
is replaced by x°(r, i"';w), Eq. (||) yields the exact ex- 
change energy per particle at point r, Sx [n] (r) . The RPA 
exc[n(v) is derived by— simply setting the xc kernel of 
Eq. (1^) equal to zero.Ea In the so-called time-dependent 
local-density approximation (TDLDA)l3 or, equivalently, 
adiabatic local-density approximation ( ALDA) , the exact 
long-wavelength limit of the static {uj = 0) xc kernel is 
used, i.e.. 



J X 



\,ALDA 



[n](r,r';w) 



dn{T) 



5{r 



(14) 



w^^""*-^[n(r)] being the xc potential of a uniform electron 
gas of density ri(r). 

For a uniform electron gas of density n, introduction 
of Eq. (0) into Eq. (0) and then Eq. ^ into Eq. (|) 
yields 



dq 
(2^ 
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dujx^{q,i^) - 1 
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where v{q) and x^{1t^) represent three-dimensional 
Fourier transforms of the Coulomb potential and the 
density-response function, respectively. If the interact- 
ing density-response functicip-|X^(g;^) replaced by the 
Lindhard function x^{lT'-^)f^ Eq. (15) yields the exact 



exchange energy per particle of a uniform electron gas. 



and the LDA exchange surface energy 



LDA 



3e^ 



^•^n dz' p{z') l-pi/3(z') 



(16) 



(17) 



where p{z) — n{z)/n represents the normalized electron 
density, and z' — 2qp z. 
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For a jellium slab, which is translationally invariant 
only in a plane perpendicular to the z axis, one finds 
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dz' z'; q||) 
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where v{z,z';q^^) and x'*'(z, z'; , w) represent two- 
dimensional Fourier transforms of the Coulomb poten- 
tial and the density-response function, and q|| is a wave 
vector parallel to the surface. 

To find x^i^: ^'i Q\\ i we follow the method described 
in Ref. We first assume that|_a(z) vanishes at a dis- 
tance zo from either jellium edgejH3 and expand the wave 
functions (z) in a Fourier sine series. We then introduce 
a double-cosine Fourier representation for the density- 
response function, which allows us to transform Eq. jl^ ) 
into a matrix equation that is solved numerically. The 
integrals of Eqs. (||) and ( p^ ) over the coordinate normal 
to the surface can be performed analytically, and we find 
an explicit expression for the xc surface energy, in terms 
of the eigenvalues ei and the Fourier coefficients of the 
eigenf unctions (l)i{z) (see Appendix A). 



III. THE INFINITE BARRIER MODEL 

In this section we consider the simplest possible mi- 
croscopic model of the jellium-slab, namely the so-called 
infinite barrier model (IBM) .L^JThis non-selfconsistent 
model has been widely discussed, EiTEj because of the sim- 
plicity of the one-electron wave functions. Nevertheless, 
we are not aware of any theoretical description of the ex- 
trapolation, within this model, of finite-slab calculations 
to the infinite- width limit. Since this issue may be im- 
portant, in view of its potential application to the case 
of more complex self-consistent calculations, we present 
below a detailed analysis of QSE on the surface energy 
of an electron system confined in an infinite well. 

In the IBM, the Kohn-Sham effective one-electron po- 
tential Veff[n]{z) is replaced by infinitely high potential 
walls at a distance zq outside the jellium slab, 

C 0, -a/2 - Zo < a/2 + zq 
Veffiz) - <^ (19) 
I oo, elsewhere, 

with Zo chosen so as to ensure charge neutrality, i.e., 

/•a/2+ZQ 

po{a + 2za) = 2 dz p{z) = a, (20) 



Po representing the normalized electron density averaged 
over the region between the two infinite barriers. 

Hence, the normalized Kohn-Sham one-electron wave 
functions and energies are 



In 



and 
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with d = a + 2zq. Introduction of these eigenfunctions 
and eigenvalues into Eq. (^ yields the IBM electron 
density, and no attempt is then made to alter Veff[n\{z) 
using this computed density. 

In the infinite- width (iw)|Jimit and for z > 0, Eq. (|^) 
yields the well-known resultu 

[piz)]^n, = [l + 35-3(z cosz-sinz)] e(-z), (23) 

where z — 2qp{z — d/2). By carrying out the integration 
of Eq. (20) one then easily finds 



and 



[^o].„ = (3/16)A, 



(24) 
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where the dimensionless parameter x — 2d/\p has been 
introduced, \p = 2'iT/qF being the Fermi wavelength. 
Instead, for a finite slab one finds 
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and 
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with the highest occupied level Im being the largest in- 
teger less than or equal to x. 

The dependence of the average electron density po on 
the parameter obtained from Eq. (p6|), is shown in 
Fig. 1. This figure clearly shows the oscillatory character 
of this quantity, which is a consequence of the so-called 
QSE first reported by Schulte.E3 As d increases, new sub- 
bands for the z motion become occupied. One easily sees 
that such a new subband falls below the Fermi level every 
time d is a multiple of Ai?/2, i.e., every time x is integer. 
When a new subband is pulled below the Fermi level, the 
parallel Fermi sea built upon the newly occupied subband 
acquires more electrons, thereby increasing p^. However, 
this effect is eventually overcome by the fact that all the 
subbands for the z motion get deeper with increasing film 
thickness. When x has increased by unity, a new subband 
begins to be filled and a new oscillation begins. 

We define the threshold value [po]i of the normalized 
average electron density po corresponding to a jellium 
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slab with x integer. By using the Euler-Maclaurin sum- 
mation formula,Ej Eq. ( p6|) with x = n yields a local 
minimum, 



3/1 



(28) 



We also define [po]2 corresponding to a jellium slab with 
X = n + 1/2, and find 



N2 = l-7^(l-i 



4x 
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(29) 



As X increases, both Eqs. (|28|) and ( |29| ) approach the 
infinite- width asymptotic behavior of Eq. (p^. Although 
[po]i and [pq\2 represent the actual value of po only at a; = 
n and x = 7T.-I-1/2, respectively, they are exhibited in Fig. 
1 for all values of x. This figure shows that Eqs. (|8|) and 
( p9| ) represent enveloping lines of the actual oscillating 
behavior of po- Hence, the amplitude of the oscillations 
is 



[Po\2 - [poll 
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(30) 



thus decaying linearly with the thickness of the slab. 
Moreover, the infinite- width limit of Eq. ([24|) is obtained 
from the enveloping lines of Eqs. (|2^) and (|29|), as follows 
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[poll + 2 [po] 
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Fig. 2 shows that a similar oscillatory behavior is ex- 
hibited by the normalized electron density p{z) at the 
center of the slab (2 = 0): 



21 - 1 



(32) 



where Z^.j is the largest integer less than or equal to [x - 
l)/2. Also shown in Fig. 2 are the enveloping lines 

1 



[p(0)]i = 1 - - 



and 



1 
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(33) 



(34) 



which represent the actual value of p{Q) for x odd and 
even, respectively. The infinite- width limit [/c(0)]j^ — 1 
is obtained as in Eq. (|3T|), with p(Q) instead of po- 

At the threshold width for which the nth subband for 
the z motion is first occupied {x = n), the normalized 
average electron density po is smaller than in the infinite- 
width limit, i.e., [po]i < [Po]iu,, which requires that the 
positive background be cut-off further from the infinite 
wells, i.e., zq is a local maximum. This is illustrated in 
Fig. 3, where the dependence of the ratio zq/ [zq]^^ on 



the parameter x is exhibited, as obtained from Eq. 
Also plotted in this figure are the enveloping lines 



and 



[2:0] 2 = [zo]^w 
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(36) 



which coincide with the actual value of zq for x integer 
{x — n) and for x = n + 1/2, respectively. An inspection 
of Eqs. ( ^5| ) and ( |3^ ) easily shows that the infinite- width 
limit of Eq. (^5|) is obtained as in Eq. (^), with po 
replaced by zq. 



A. Kinetic surface energy 

An explicit expression for the surface-energy term Us 
is obtained by substituting Eqs. ( p^ and ( p^ into Eq. 
(i). One finds 



= K],^ 10 £ 
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where 
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is the often-derived value of the IBM surface-energy term 
(Ts of a semi-infinite mcdium.cl For x integer, the Euler- 
Maclaurin summation formula now yields 



1 
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(-i)l 


3a; 





while for X = n + 1 /2 one finds 
1 



'sJ2 



The dependence of the ratio a 



(39) 



on the param- 
7), is exhibited in Fig. 




eter x, as obtained from Eq. 
4. Also plotted in this figure are [(Js]i and [as]2- While 
[(Ts]i represents the actual surface-energy term CTs at the 
local minima {x = n), local maxima slightly deviate from 
X — n + 1/2, especially for the smallest values of x. For 
large values of a;, Eq. ( pO[ ) approaches the upper envelop- 
ing line of the true oscillating behavior of as- Further- 
more, one finds 



[as],+2[as], 



+ 0(x-3), 



(41) 



where the first term quickly approaches the actual 
infinite- width limit of Eq. (pq), as shown in Fig. 4. 
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B. Exchange surface energy 

The LDA exchange surface energy is obtained by sim- 
ply introducing the normaUzed electron density into Eq, 
(p^). In the infinite- width limit p{z) is that of Eq. (2 



[a.,],^ = 4.0700 X 10- 



(45) 



and numerical integration leads to the following resultl 
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6.3179 X 10 
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(42) 



as reported in Ref. ^. 

For the actual IBM electron density of a finite slab of 
thickness d — xAf/2, introduction of Eq. (|^) into Eq. 
( p7| ) yields the ratio cr^^"*/ [o'i'^'*] shown in Fig. 5, as 
a function of x. As in the case of the kinetic surface en- 
ergy, we observe a considerable quantum-size effect, even 
for slab widths as large as five times the Fermi wave- 
length. The amplitude of the oscillations is also found to 
decay linearly with the thickness of the slab. However, 
explicit expressions for the enveloping lines [c^^"^] ^ and 
[fj^^'^] 2 are not available; hence, for each integer value 



which is in exactj-agreement with the result obtained by 
Harris and JonestI for the semi-infinite medium. 

The extrapolated exact exchange surface energies, as 
obtained from Eq. (44) with a replaced by a^, are rep- 
resented in Fig. 6 versus x (only integer values of x have 
been considered) by open squares, showing that it ap- 
proaches the infinite- width limit of Eq. (^). Although 
the extrapolation procedure dictated by Eq. ( ^ ) does 
not converge as rapidly for ax as for a^^^, the numerical 
error introduced by the use of this procedure to evaluate 
ax is still found to be within 0.7% for slab widths of only 
5 times the Fermi wavelength {x = 10). 



^LDA] 



of X {x — n) we take [a^ j ^ 
mate [ct^^'*] ^ as follows 
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where [a^^^]^ = a^^^{n - 1/2) and [a^''^ = 
a^^^{n + 1/2). Assuming that the infinite- width limit 
can be extrapolated as in Eq. (^l|), we use the relation 



(44) 



with a replaced by a^^^. 

This procedure yields the results represented in Fig. 
5 by open circles, which quickly approach the actual 
infinite- width limit of Eq. (|4^). The normalized LDA 
exchange surface energies obtained from Eq. (^) are 
also shown in Fig. 6 (open circles), but now for val- 
ues of x up to 100 and together with local maxima and 

minima, [-.^^^] ^ / k^^^] and k^^^] / [a,^^^] 
respectively (dashed lines). While the amplitude of the 
quantum-size effect amounts, for slab widths as large as 
50 times the Fermi wavelength {x — 100), to 1% of the 
actual infinite- width limit, for slab widths of only 5 times 
the Fermi wavelength {x ~ 10) the numerical error intro- 
duced by the use of Eq. (^) is found to be within 0.02%. 

Within the IBM, we have also calculated exact ex- 
change surface energies by introducing Eqs. (|l^) and 
(|l|) into Eq. d) with x°(g,w) and x°{z,z';qii,u;) in- 
stead of x'*'(9,w) and x'*'(z, z'; (?|| , w). We have found an 
oscillatory behavior of this exact contribution to the sur- 
face energy that is similar to that shown in Fig. 5, and 
have tested the approximation of Eq. (|^ (with a re- 
placed by ax) for slab widths as large as 100 times the 
Fermi wavelength {x 200). We have found 



C. Exchange-correlation surface energy 

For a given electron density, the LDA xc surface energy 
is obtained by introducing into Eq. (H) the xc energy per 
n) and approxi- particle of a uniform electron gas, e™*-^, which we have 
calculated, with use of either the Wigner intecaolation 
formulaJ23 the Perdew-Zunger parametrizatiorO of the 
QMC calculation of Ceperley and Alder,E3 or Eq. ( p^ ) 
with the RPA interacting density-response function. We 
have found quantum-size effects and the extrapolation 
of finite-slab calculations to the infinite-width limit to be 
similar to those exhibited by kinetic and exchange contri- 
butions to the surface energy. Hence, for each value of 
and a given number n of occupied subbands, we have con- 
sidered three different values of the parameter x: n — 1/2, 
n, and n-|-l/2, and have used Eq. (|j) with a replaced by 
a^^^ . When either the Wigner interpolation formula or 
the RPA are used for e"™-'^, the extrapolation procedure 
dictated by Eq. ( ^ ) yields LDA xc surface energies that 
are in agreement with those reported for a semi-infinite 
medium in Refs. |^ and ^, respectively. Since the Wigner 
correlation energy changes very little with the electron 
density, as compared with RPA or Perdew-Zunger corre- 
lation energies, it yields LDA xc surface energies that are 
too small. 

For each and given values of a;, we have also com- 
puted non-local xc surface energies from Eq. (|5|), as ob- 
tained with the xc energy of Eqs. (|l^) and ( |18D com- 
puted in the RPA. We have found that the extrapo- 
lation of these finite-slab calculations to the infinite- 
width limit is accurately described by Eq. (0), but 
now with a replaced by axe- Table I shows our extrap- 
olated RPA correlation-only and exchange-correlation 
surface energies, ac and axe, together with the RPA- 
based HM^ surface energies, a^^^ and a^^^, described 
above .l3ej While the local-density approximation over- 
estimates, within the IBM, the exact exchange surface 
energy by 55% [see Eqs. (^ and (^5|)] for all values of 
Ts , it largely underestimates the RPA correlation surface 
energy by 78% at = 2 and 64% at Ts = 6. Large and 
opposite separate corrections to the LDA for exchange 
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and correlation nearly compensate, as discussed in Refs. 
^ and ^, and yield LDA xc surface energies that are, over 
the entire metallic density range, within about 6.5% of 
their non-local counterparts. 

For comparison. Table I also shows the extrapolated 
LDA correlation-only and exchange-correlation surface 
energies, c^^"^ and a^^^ , that we have obtained with 
use of the Perdew-Zunger parametrization for e"™^ . Al- 
though RPA and Perdew-Zunger correlation energies of 
a uniform electron gas differ, for ~ 0-10, by ~ 0.6- 
0.3 eV [the second-order exchange energy not present in 
the RPA is known to yield ~ 0.6 eV in the high-density 
limit {r, 0)], the quantity {e™^/[n(0)] - en- 
tering Eq. depends little on whether the RPA or the 
Perdew-Zunger parametrization is used, thereby yielding 
LDA xc surface energies that are comparable. 

IV. SELF-CONSISTENT CALCULATIONS 

Due to the simple nature of the IBM one-electron wave 
functions, we have used this model as a proving ground 
for the calculation of both local (LDA) and non-local xc 
surface energies, and also for devising the extrapolation 
procedure dictated by Eq. (|^. Nevertheless, the IBM 
yields unrealistic xc surface energies that are too small, 
especially at the lowest values of rg that we have explored, 
since at these high densities the infinite barrier constrains 
the electrons too close to the jellium background. 

In this section, we apply Eq. ( ^4|) in the framework 
of a self-consistent microscopic model of the jellium slab. 
While in the IBM new subbands for the z motion become 
occupied every time d is a multiple of Af/2 {x integer), 
for an arbitrary effective one-electron potential new sub- 
bands are expected to enter as 2a/ Xp (not necessarily 
integer) is increased by approximately one. Accordingly, 
surface energies are still oscillatory functions of the slab 
thickness a, the period of the oscillations being approx- 
imately equal to Xp/^, and the amplitude decaying ap- 
proximately linearly with a. 

A. LDA orbitals 

Self-consistent calculations of jellium surface energies 
of a semi-jnfinite medium were carried out by Lang 
and KohnB with use of the LDA for exchange and cor- 
relation. They replaced the actual VxdnKz) entering 
Eq. (^ by the LDA xc potential and computed the 
LDA-^c surface energy from Eq. (||) using the Wigner 
forniE3 for the correlation energy per particle of the 
uniform electron gas, e" ™-^. We have tested the ap- 
proximation of Eq. m4) (with a replaced by our self- 
consistent a^^^), for slab widths as large as ~ 10 times 
the Fermi wavelength, and for = 2 we have found 
lai:?"^] . = 3256 erg/cm^, in agreement with the calcu- 
lation reported by Lang and Kohn.EJ Our self-consistent 



calculations of a^^^ / [i^xc^]i^ are plotted in Fig. 7, 
together with the enveloping lines ^ / [(^xF'^] iw 

and [cr^c^^^] 2 / [o'xi''^] , as a function of the slab width 
(in units of Af/2). Also represented in this figure are the 
extrapolated LDA xc surface energies (open circles), as 
obtained from Eq. ( |4^ ) for those threshold slab widths a„ 
for which the nth subband is first occupied. One clearly 
sees that self-consistent density profiles yield exactly the 
same oscillating behavior as obtained within the IBM, 
but now QSE are smaller and the extrapolation formula 
of Eq. (^) yields a result that for a slab width as little 
as ^ 3.5 times the Fermi wavelength [n = 7j)-is within 
0.01% of the actual infinite-width calculation.^ 

Now we focus on our RPA-based local and non-local 
calculations of the surface energy. As in the calculations 
of Lang and Kohn, we replace the actual Vxc[n\{z) en- 
tering Eq. (|^) by the LDA xc potentialjlj but with the 
Perdew-Zunger parametrization for the xc potential of 
the uniform electron gas. Local and non-local calcula- 
tions of the xc surface energy are then carried out within 
one and the same density-response framework, i.e., by ei- 
ther introducing Eq. ([l5| ) into Eq. (@) or by introducing 
both Eqs. ( p^ ) and ( [l^ into Eq. (P), with the use of 
the interacting RPA density-response function. As this 
framework is devoid of any ambiguities in the treatment 
of the many-body problem, it clearly serves as a measure 
of the actual error introduced by the LDA. 

In Table II, we show the infinite-width limit of the 
various contributions to our RPA-based local (LDA) and 
non-local jellium surface energies, as obtained from Eq. 
(Q) for slabs with n = 12 (a ~ 5 — 6 A_f).e3 As in the case 
of IBM orbitals, self-consistent calculations show that the 
LDA systematically overestimates/underestimates the 
exchange/correlation surface energies, while the local xc 
surface energies differ little from their nonlocal counter- 
parts. Table II shows that the percentage error between 
local and non-local exchange surface energies monotoni- 
cally decreases, within the present model, from a mini- 
mum of 16% for the slow varying density at = 2 to 
a maximum of 95% for the more rapidly varying density 
at Ts = 6. The percentage error between our local and 
non-local xc surface energies is found to vary from 1.8% 
at rs = 2 to 3.4% at — 6. 

Self-consistent electron densities vary slowly as com- 
pared to the IBM, especially in the case of small values 
of r^, as shown in Fig. 8. As a consequence, we find: 
(a) quantum-size effects are smaller [for slabs with 12 
occupied levels the amplitude of the oscillations in both 
a^^^ and axe ranges from ~ 1% at rs = 2 to ~ 5% at 

=6]; (b) all contributions to the surface energy are 
larger; and (c) the LDA is significantly better, over the 
entire metallic range of densities (r^ ~ 2 — 6), than for 
the IBM electron-density profile. 

The impact on the surface energy of the short-range 
correlations built into the many-body kernel of Eq. (|l^), 
not present in the RPA, was investigated in Ref. |45| 
within the ALDA. While introduction of the ALDA ker- 
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nel slightly decreases the LDA xc surface energy [the 
ALDA-based LDA surface energies agree closely with the 
LDA quantum Monte Carlo calculations of Acioli and 
Ceperleya], ALDA short-range corrections to the non- 
local RPA surface energy are found to be positive. As a 
result, the presence of the xc kernel enhances, within the 
ALDA, the discrepancy between local and non-local xc 
surface energies, but the error introduced by the LDA is 
still, within ALDA, of the order of 50% smaller than the 
error one would impute to the local apppaximation on 
the basis of the non- local FHNCA results,E3 particularly 
in the high-density region (r<, = 2). 

B. Hartree and exchange-only orbitals 

RPA-based local and non-local surface energies, as ob- 
tained with the use of Hartree orbitals, i.e., with the xc 
potential W2:c[f^] entering Eq. set equal to zero, were 
reported in Ref. Hartree orbitals are more delocalized 
than the more realistic LDA orbitals (with Vxc[n']{z) = 
the work function is far too small and even negative 
at low values of r^) and they yield a too slow varying 
electron-density profile, as shown in Fig. 8. This leads 
to surface energies that are too large, relative to those 
obtained with the use of LDA orbitals. We also find that 
the impact of xc on the one-electron states is much more 
pronounced than the impact of the short-range correla- 
tions built into the many-body kernel of Eq. (^3|) . 

For comparison, we have also computed kinetic, elec- 
trostatic and exchange surface energies with an LDA 
exchange-only effective one-electron potential Vx[n\{z), 
and have obtained the infinite- width extrapolated results 
shown in Table III. A comparison of the numbers pre- 
sented in Tables H and III shows that the sum of kinetic, 
electrostatic and exchange surface energies is nearly in- 
sensitive to whether an LDA exchange-only or xc effec- 
tive potential is used in the evaluation of the one-electron 
Kohn-Sham orbitals. Moreover, this sum (labeled (Jhf 
in Table HI) is expected to be jffiry close to the ex- 
act Hartree-Fock surface energies.O Our exchange-only 
surface energy is for rg = 2.07 in excellent agreement 
with the Hartree-Fock surface cncpgy of —1273 erg/cm^ 
reported by Krotscheck and Kohn.E3 Also, our exchange- 
only surface energies fall for > 3 belo^^^he variational 
upper bound reported by Sahni and Mao 

V. SUMMARY AND CONCLUSIONS 

We have presented self-consistent-field calculations of 
the surface energy of a bounded electron gas. First of 
all, we have considered the infinite-barrier model of the 
surface. Within this model we have discussed the origin 
of quantum-size effects, and have devised a formula that 
accurately extrapolates our finite-slab calculations to the 
infinite- width limit. This formula may be useful to the 



case of more complex slab calculations involving the use 
of the band structure and/or quantum Monte Carlo sim- 
ulations. 

Subsequently, we have reported full self-consistent non- 
local calculations of the surface energy, by employing ei- 
ther LDA or the less realistic Hartree orbitals to con- 
struct the RPA density-response function. The unam- 
biguous nature of the comparison of local versus non-local 
surface energies made possible by our self-consistent cal- 
culations leads us to the conclusion that the local-density 
approximation is accurate, within the RPA treatment of 
correlation. 

Our non-local RPA surface energies, as computed 
with the use of realistic LDA orbitals, are in reason- 
able agreemottt with those obtained by Zhang, Langreth, 
and Perdfiy^ using the non-local Langreth-Mehl xc 
functional and also with those obtained within the re- 
cently developed jneta-generalized-gradient approxima- 
tion (meta-GGA),E3 which makes use not only of the local 
density and its gradient but also p£ the orbital kinetic- 
energy density. Moreover, a GGAcj for the short-range 
correlation not present in the RPA has yielded, a small 
and negative contribution to the surface energy ,tZI thereby 
suggesting that the RPA may be a much better approxi- 
mation for the changes in correlation energy upon surface 
formation than for the total correlation energy. A recent 
wave-vector, interpolation as a long-range correction to 
the GGA,Ej which has predicted surface energies that 
are smaller but close to our non-local RPA results, has 
indicated that the LDA should be accurate even within 
an exact treatment of electron correlation. This is in con- 
trast with non-local correlation surface energies reported 
in Refs. ^ and ^ which arc significantly higher than our 
non-local RPA correlation energies and the local-density 
approximation. 
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APPENDIX A: 

Here we give an explicit expression for the xc surface 
energy of a jellium slab, in terms of the eigenvalues Si 
and the Fourier coefficients of the eigenfunctions 4'i{z). 
It is also demonstrated that there is no contribution to 
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the xc surface energy coming from large values of the mo- 
mentum transfer parallel to the surface, so the xc energy 
per particle at point z is devoid of divergences. 

We first introduce the following double-cosine Fourier 
representations for the density-response function and the 
Dirac (5 function: 



X(2,z';g||,w) = ^ ^Y^Xmn{q\\,'jJ) cos(m7rz) cos(n7rz') 

(Al) 



m— n— 



and 



b{z — z') = cos(ri7rz) cos(r7,7rz'), (A2) 



ri=0 



where d = a + 2zo, z = z/d+ 1/2, and /Xm is defined by 
the equation 



1 for m = 0, 



2 for TO > 1. 



(A3) 



The electron density can also be expressed in a cosine 
representation 



i{z) — Um cos(TO7rz). 



(A4) 



m— 



Eq. m yields 



(A5) 



1=1 



where 



G™' = d / cos(TO7rz) 0/(z) (/);/ (z), (A6) 
Jo 

or, in terms of the Fourier coefficients bis of (f>i{z), 

^ OQ OO 

s=l s' = l 

(A7) 



Then, substituting Eq. ^ into Eq. (|), we obtain 
the following expression for the xc surface energy: 



crxc = ^ ^ ^ ""m(9||) -~naea,c{n), (A8) 



oo oo 



m— n— 



where 

0^mn{q\\) = "^mn (g||) 



"'0 



- dX dwxmnin^i^) 



1=1 



i'=i 



and 

Vmn{q\\) 



27re2 



g2 (TO7r/d)2 

-[i + (-ir 



2d 



,„i 911 [i-(-ir 



(Q| + (n7r/d)2 



(AlO) 



As for the Fourier coefficients Xmn(9|h^); Eq. ( [l^ ) 
yields 

X |Aw,„/„'(g||) -f/^';;f,[n](q||,w)|x^„(g||,w), (Al 



1) 



where the Coulomb-interaction Fourier coefficients 
Vmni,q\\) are given by Eq. (AlO). An explicit expres- 



sion for the coefficients Xmni^h^) can be found in Ref. 
H- fmn^[^]{'l\\j^) are the Fourier coefficients of the xc 
kernel. In the RPA these coefficients are simply zero. In 
the ALDA, 



j-xc,X,ALDAi 
J ran V 



d r d~z '^''^^^''^ 



dn 



n=n{z) 



X cos(TO,7rz) cos(n7rz). (A12) 



If the Coulomb correlations entering Eq. (All) are 
ignored altogether (A = 0), the integrations over A and uj 
in Eq. ( |A9| ) can be performed analytically to find, after 
some algebra, the following result: 



Im oo 

EE 



1=1 i'=i 

[Fu4.q\\)-{EF-ei)]G';^,GTi,, (A13) 



where 

Fii'{q\\) = ei sgn[a,;/(g||)] +Q[bw{q\\)] 
a'iv{q\\)y'hiv{q\\) 2 



irq 



■ ei arctan 



aw 



'^n'{q\\) = ^ql-\{si-e,l 



and 



biv{q\\) = 2(?jjei - afi, 
In the limit as q\\ — s- cx), we find 



9|| 



lim Fiv{q{) = {Ep - £;)■ 



(A14) 



(A15) 



(A16) 



(A17) 



(A9) 



As in this limit the interacting density-response func- 
tion approaches its non-interacting counterpart, this re- 
sult shows, explicitly, that there is no contribution to the 
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xc surface energy coming from high values of the paral- 
lel component of the momentum transfer, and that the 
singular term S{z — z') in Eq. ( |l^ ) gives no divergences. 
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FIG. 1. The solid line represents the normalized average 
electron density po, as obtained from Eq. (|2^), versus the pa- 
rameter X = 2d/\F- Dashed lines represent [po]i and [po]2i as 
obtained from Eqs. ( p^ ) and (p9|). The dotted line represents 
the infinite- width limit of Eq. (pj) . 

FIG. 2. As in Fig. 1, for the normalized electron den- 
sity at the center of the slab, obtained from Eqs. (|32[), (p3[), 
and (^^. The dotted line represents the infinite- width limit, 
[P(0)],„ = 1. 



TABLE I. Non-local RPA correlation (ctc) and ex- 
change-correlation ((Tj^c) surface energies obtained in the infi- 
nite-barrier model of the surface, and their local counterparts. 
Also shown are the LDA xc surface energies {o-^^Jr'^ ) obtained 
by using the Perdew-Zunger interpolation formula for the cor- 
relation energy of a uniform electron gas. Units are erg/cm^. 
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LDA' 
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LDA' 
^ xc 
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1357 


2.07 


617 
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99 


31 


26 


198 


185 


180 
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53 


18 


15 


104 


97 


94 


6.0 


33 


12 


10 


62 


58 


56 



TABLE II. Kinetic ((Ts), electrostatic (cfcs), exchange 
((Tj;), non-local RPA xc ((Txc), and total (ct) surface ener- 
gies and their local counterparts, as obtained with the use 
of Perdew-Zunger parametrized LDA orbitals. Units are 
erg/cm^. 
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TABLE III, 


, As 


in Table II, but with use of exchanj 


^e-only 


LDA orbitals. 
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32 
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FIG. 3. As in Fig. 1, for zq/ XMi^ obtained from Eqs. (|^, 
(H), and 

FIG. 4. As in Fig. 1, for Osj [fsJi^ obtained from Eqs. 
(0), (!§), and (§). 

FIG. 5. The solid line represents the normalized LDA ex- 
change surface energy, cr^^ '^/ [cr^^^] , , as obtained from 
Eq. ( p^ for the IBM electron density of a finite slab, versus 
the parameter x — 2d/\F- Dashed lines represent the en- 
velopes and / [a^^] . The 

extrapolated normalized LDA surface energies, as obtained 
from Eq. (^) with a replaced by a^^^, are represented by 
open circles. 

FIG. 6. The solid lines 

represent the envelopes [ux]^ / ['^x]^^ and [cr^^'^jj / k^li^, of 
the normalized exact exchange surface energies, which have 
been obtained, within the IBM, by introducing Eqs. (^^ and 
( p^ into Eq. (^ with ^{q,u>) and -x^\z,z' \q^^,u) instead 
ol X'^{q,Lo) and x^i^^ ^' The corresponding LDA en- 
veloping lines (also shown in Fig. 5) are represented by dashed 
lines. Open squares and circles represent extrapolated exact 
(squares) and LDA (circles) exchange surface energies, as ob- 
tained from the procedure dictated by Eq. (^. 

FIG. 7. The solid line represents the normalized LDA xc 
surface energy, o^^^ j [cr^^'*] . , as obtained from Eq. (^) 
for the Lang-Kohn electron-density profile of a finite slab, ver- 
sus the slab width (in units of Af/2). Dashed lines represent 
the envelopes {a^^°^\ ^ / {a^^°^\ and [a,^^^] ^± / {a^^°^\ 
The extrapolated normalized LDA surface energies, as ob- 
tained from Eq. ( ^ ) with a replaced by o'^i''*, are repre- 
sented by open circles. 
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FIG. 8. The normalized electron-density profile of a jellium 
surface with Vg = 2.07, as obtained with use of Perdew-Zunger 
parametrized LDA orbitals (solid line) , Hartroo orbitals (dot- 
ted line), and within the IBM (dashed line). We note that 
the normalized IBM electron-density profile, with the z coor- 
dinate measured in units of the Fermi wavelength, is indepen- 
dent of Tg. In this figure, the jellium edge has been taken at 
z = G. 
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